Scipy

发表于 | 分类于 ， | | 本文总阅读数 次
| 字数统计: 2.1k | 阅读时长 ≈ 10
Scipy

0

Scipy介绍

  Scipy是一个用于数学、科学、工程领域的常用软件包，可以处理插值、积分、优化、常微分方程数值解的求解、图像处理、信号处理等问题。它用于有效计算Numpy矩阵，使Numpy和Scipy协同工作，高效解决问题。

Scipy特点

  Scipy支持大多数工程数学运算。
  Scipy每一个子模块都可以完成一类功能。
  Scipy中的函数类似于MATLAB中的函数，使用方便。

Scipy应用

Scipy常数模块

1
2
3
4
5
6
7
8
9
10
11
12
13
from scipy import constants as C

# C.c 光速常数
C.c

# C.h 普朗克常数
C.h

# C.mile 英里
C.mile

# C.pi 圆周率π
C.pi

1

Scipy特殊函数模块

gamma，gammaln方法

1
2
3
4
5
6
7
from scipy import special as S

# S.gamma(n) 计算Γ(n)的值
S.gamma(4)

# S.gammaln(n) 计算ln|Γ(n)|的值，避免Γ(n)过大
S.gammaln(4)

2

log1p方法

1
2
3
4
5
6
from scipy import special as S
import numpy as np

# S.log1p(n) 计算ln(n+1)的值，使其可以计算很小的数
S.log1p(np.e - 1)
S.log1p(1e-10)

3

Scipy拟合优化模块

fsolve方法

1
2
3
4
5
6
7
8
9
from scipy import optimize as O
import math

def f(x):
    x0, x1, x2 = x.tolist()
    return [5 * x1 + 3, 4 * x0 ** 2 - 2 * math.sin(x1 * x2), x1 * x2 - 1.5]

# O.fsolve(f, init) 求非线性方程组的解，f为方程函数，init为初始迭代值
result = O.fsolve(f, [1, 1, 1])

$$ \begin{cases} 5x_1 + 3 & = 0 \\ 4{x_0}^2 -2\sin(x_1 \cdot x_2) & = 0 \\ x_1 \cdot x_2 -1.5 &= 0 \end{cases} $$
4

leastsq方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from scipy import optimize as O
import numpy as np
import matplotlib.pyplot as plt

x = np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
y = np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])

def residuals(p):
    k, b = p
    return y - (k * x + b)

r = O.leastsq(residuals, [1,0])

k, b = r[0]
y_new = x * k + b
print('k=', k, 'b=', b)

plt.scatter(x, y)
plt.plot(x, y_new)
plt.show()

5

Scipy线性代数模块

solve方法

1
2
3
4
5
6
7
8
from scipy import linalg as L
import numpy as np

a = np.array([[2, 2, -1], [1, -2, 4], [5, 8, -1]])
b = np.array([[6], [3], [27]])

# L.solve(A, b) 求线性方程组Ax = b的解
x = L.solve(a, b)

6

eig，svd方法

1
2
3
4
5
6
7
8
9
10
from scipy import linalg as L
import numpy as np

a = np.array([[-2, 1, 1], [0, 2, 0], [-4, 1, 3]])

# L.eig(array) 求array的特征值和特征向量
m, x = L.eig(a)

# L.svd(array) 求array的奇异值分解
u, sigma, v = L.svd(a)

7

Scipy统计模块

norm类，stats，rvs方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from scipy import stats as ST

# ST.norm(loc=0, scale=1) 获取偏移为loc(默认为0)，标准差为scale(默认为1)的正态分布(还可以定义其他的分布)
norm_ = ST.norm(loc=1, scale=2)

# obj.stats() 获取obj分布的期望和方差
mean_norm, var_norm = norm_.stats()

# obj.rvs(size=shape) 获取大小为shape的obj分布的随机抽样
x = norm_.rvs(size=(100, 100))
mean_x = np.mean(x)
std_x = np.std(x)

# obj.pdf(x) 获取x处的概率密度函数
pdf_1 = norm_.pdf(1)

# obj.cdf(x) 获取x处的分布函数
cdf_1 = norm_.cdf(1)

8

rv_discrete类，stats，rvs方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from scipy import stats as ST

# ST.rv_discrete(values=(x, p)) 自定义离散概率分布，x为可能的取值，p为对应的概率
discrete_ = ST.rv_discrete(values=([1, 2, 3, 4, 5, 6], [0.75, 0.05, 0.05, 0.05, 0.05, 0.05]))

discrete_.stats()

x = discrete_.rvs(size=(100, 100))
mean_x = np.mean(x)
var_x = np.var(x)

plt.plot(x, y)

# plt.annotate(formulation, xy, xytext, fontsize, arrowprops) 在xy附近添加注解formulation，xytext为注解的位置，fontsize为注解的大小，arrowprops为箭头类型
plt.annotate(r'$x^2-0.5$', (0, -0.5), xytext=(+0.25, 0), arrowprops=dict(arrowstyle='->', connectionstyle='arc3, rad=0.2'))

plt.show()

9

binom.pmf方法

1
2
3
4
from scipy import stats as ST

# ST.binom.pmf(list, n, p) 进行n次二项分布实验，出现的概率为p，计算出现list中对应值的概率
x = ST.binom.pmf([0,1,2,3,4], 3, 0.8)

10

Scipy积分模块

quad方法

1
2
3
4
5
6
from scipy import integrate as I

I.quad(f, min_lim, max_lim) 计算函数f的积分(f可以是自定义函数也可以为lambda表达式)，并返回计算产生的误差，积分下限为min_lim，积分上限为max_lim

pi_half, err = I.quad(lambda x:(1 - x ** 2) ** 0.5, -1, 1)
print('pi:', pi_half * 2, 'err:', err)

$$\int_{-1}^{1} \sqrt{1-x^2}, dx = \frac{\pi}{2}$$
11

dblquad，tplquad方法

1
2
3
4
5
6
7
8
9
10
from scipy import integrate as I

f = lambda x, y : (1 - x ** 2 - y ** 2) ** 0.5
f_y = lambda x : (1 - x ** 2) ** 0.5

# I.dblquad(f, x_min_lim, x_max_lim, y_min_lim, y_max_lim) 计算函数f的二重积分，并返回计算产生的误差，x积分下限为x_min_lim，积分上限为x_max_lim，y积分下限为y_min_lim，积分上限为y_max_lim
res, err = I.dblquad(lambda x, y : (1 - x ** 2 - y ** 2) ** 0.5, -1, 1, lambda x : -1 * f_y(x), f_y)
print('res:', res * 1.5, 'err:', err)

# I.tplquad(f, x_min_lim, x_max_lim, y_min_lim, y_max_lim, z_min_lim, z_max_lim) 计算函数f的三重积分，并返回计算产生的误差，x积分下限为x_min_lim，积分上限为x_max_lim，y积分下限为y_min_lim，积分上限为y_max_lim，z积分下限为z_min_lim，积分上限为z_max_lim

$$\int_{-1}^{1} \int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}} \sqrt{1-x^2-y^2}, dydx = \frac{2}{3}\pi$$
12

Scipy插值模块

interp1d，interp2d方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from scipy import interpolate as IP
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 2 * np.pi, 11)
y = np.sin(x)
x_new = np.linspace(0, 2 * np.pi, 51)

# IP.interp1d(x, y, kind) 对x，y进行一维插值，kind为插值函数类型，可以为'nearest'最近邻插值，'cubic'立方插值等等
f = IP.interp1d(x, y, kind='cubic')
y_new = f(x_new)

plt.subplot(121)
plt.plot(x, y, label='origin')
plt.subplot(122)
plt.plot(x_new, y_new, label='new')
plt.legend()

plt.show()

# IP.interp2d(x, y, z, kind) 对x，y，z进行二维插值，kind同interp2d

13

UnivariateSpline方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from scipy import interpolate as IP
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 2 * np.pi, 11)
y = np.sin(x)
x_new = np.linspace(-0.5 * np.pi, 2.5 * np.pi, 51)

# IP.interp1d(x, y, w, k=3, s=None) 对x，y进行一维插值，w为每个数据的权值，k为插值的阶数(默认为3)，s为曲线的平滑系数(默认为None)，和interp1d不同的是，该方法支持外推操作，即可以插值边缘点之外的部分。
y_new=ip.UnivariateSpline(x, y)(x_new)
plt.subplot(121)
plt.plot(x,y,label='origin')
plt.subplot(122)
plt.plot(x_new,y_new,label='new')

plt.show()

14

Scipy信号处理模块

medfilt方法

1
2
3
4
5
6
7
8
9
from scipy import signal as SP
import numpy as np

np.random.seed(1)

x =np.random.randint(0, 9, (5, 5))

# SP.medfilt(array, kernel_size) 对array进行中值滤波，掩模的大小为kernel_size
y = SP.medfilt(x, 3)

15

order_filter方法

1
2
3
4
5
6
7
8
9
10
11
from scipy import signal as SP
import numpy as np

np.random.seed(1)

x =np.random.randint(0, 9, (5, 5))

# SP.order_filter(array, domain, rank) 对array进行模板为domain的排序滤波，rank为第几小的值，rank=0代表最小值滤波，rank=domain.size-1代表最大值滤波
y_min = SP.order_filter(x, np.ones((5, 5)), 0)
y_mid = SP.order_filter(x, np.ones((5, 5)), 12)
y_max = SP.order_filter(x, np.ones((5, 5)), 24)

16

iirdesign，lfilter方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
from scipy import signal as SP
import numpy as np
import matplotlib.pyplot as plt

fs = 100
t = np.arange(0, 2, 1 / fs)
y = np.sin(2 * np.pi * 2 * t) + np.sin(2 * np.pi * 4 * t) + np.sin(2 * np.pi * 6 * t)

# SP.iirdesign([pass_low, pass_high], [stop_loss, stop_high], gp, gs) 设计IIR滤波器，通带频率为[pass_low × f0, pass_high × f0]，阻带频率为[0, stop_loss × f0] ∪ [stop_high × f0, ∞]，其中f0为采样频率的一半，gp为通带的最大增益衰减，gs为阻带的最小增益衰减，返回值为滤波器分子和分母的系数
b_1, a_1 = SP.iirdesign([0.05, 0.2], [0.01, 0.5], 2, 40)
b_2, a_2 = SP.iirdesign([0.1, 0.2], [0.01, 0.5], 2, 40)

# SP.lfilter(b, a, x) 计算x经过b，a滤波器的结果
out_1 = SP.lfilter(b_1, a_1, y)
out_2 = SP.lfilter(b_2, a_2, y)

plt.subplot(321)
plt.plot(t, y, label='origin')
plt.legend()

plt.subplot(322)
plt.plot(np.abs(np.fft.fft(y)), label='origin fft')
plt.xlim((0, fs / 2))
plt.legend()

plt.subplot(323)
plt.plot(t, out_1, label='out_1')
plt.legend()

plt.subplot(324)
plt.plot(np.abs(np.fft.fft(out_1)), label='out_1 fft')
plt.xlim((0, fs / 2))
plt.legend()

plt.subplot(325)
plt.plot(t, out_2, label='out_2')
plt.legend()

plt.subplot(326)
plt.plot(np.abs(np.fft.fft(out_2)), label='out_2 fft')
plt.xlim((0, fs / 2))
plt.legend()

plt.show()

17

Scipy图像处理模块

  Scipy.ndimage 是一个处理多维图像的函数库，其中又包括以下模块。
  filters 图像滤波器函数库
  fourier 傅里叶变换函数库
  interpolation 图像变换函数库
  morphology 形态学操作函数库

  图像处理有许多更强大的库，如opencv，scikit-image库，在此不做过多介绍，可以参考opencv。

Scipy小结

  通过Scipy，使用者可以仅需要几行代码，便可以完成一系列工程应用。在数据分析，实际项目中，常常需要对数据进行插值、拟合、优化，需要借助Scipy科学计算库的帮助，因此Scipy是工程研究中必不可少的帮手。

-------------本文结束感谢您的阅读-------------
0%